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ABSTRACT 

The multidimensional character of the hydrodynamics in core-collapse supernova (CCSN) cores 
is a key facilitator of explosions. A major goal over the last decade has been to elucidate how 
multidimensionality produces more favorable conditions for explosions. Unfortunately, much of this 
work has necessarily been performed assuming axisymmetry and it remains unclear whether or not this 
compromises those results. In this work, we present analyses of simplified two- and three-dimensional 
CCSN models with the goal of comparing the multidimensional hydrodynamics in setups that differ 
only in dimension. Not surprisingly, we find many differences between 2D and 3D models. While 
some differences are subtle and perhaps not crucial to understanding the explosion mechanism, others 
are quite dramatic and make interpreting 2D CCSN models problematic. In particular, we find that 
imposing axisymmetry artificially produces excess power at the largest spatial scales, power that has 
been deemed critical in the success of previous explosion models and has been attributed solely to the 
standing accretion shock instability. Nevertheless, our 3D models, which have an order of magnitude 
less power on large scales compared to 2D models, explode earlier. Since we see explosions earlier 
in 3D than in 2D, the vigorous sloshing associated with the large scale power in 2D models is either 
not critical in any dimension or the explosion mechanism operates differently in 2D and 3D. Possibly 
related to the earlier explosions in 3D, we find that about 25% of the accreted material spends more 
time in the gain region in 3D than in 2D, being exposed to more integrated heating and reaching 
higher peak entropies, an effect we associate with the differing characters of turbulence in 2D and 
3D. Finally, we discuss a simple model for the runaway growth of buoyant bubbles that is able to 
quantitatively account for the growth of the shock radius and predicts a critical luminosity relation. 
Subject headings: hydrodynamics, neutrinos, stars: interiors, supernovae: general 



1. INTRODUCTION 

Multidimensional hydrodynamic instabilities play a 
central role in the explosion mechanism of most 
core-collapse supernovae (CCSNe). Though this 
has been repeatedly demonstrated in a variety of 
contexts dHerant et all Il992t IBurrows fc Frvxelll 119931: 
Herant et al.lll994tlBurrows . et al.lll995HJanka fc Mueller! 
1996; Mar ek fc Jankal 12009). there is as yet no defini- 
tive understanding of the role of multidimensional 
effects in facilitating explosions. Important effects 
may include increased dwell times in the gain re- 
gion (|Murphy fc Burrows 2008), expansion of the shock 
due to turbulent pressure support f rom neutrino- 
driven convection dHerant et al.l 11994 IBurrows et al.l 
Il995t iMurphv et al.l I2012D and/or the standing accre- 
tion shock instability (SASI) (e.g. iScheck et ail [ 2008). 
simul taneous accretion and explosion (|Burrows et al.1 
suppression of cool ing beneath the gain region 
(|Peicha fc Thompson! [2012T ). and other still unidentified 
processes. 

With the exception of a few preliminary results in 3D 
(iBruenn et al.ll2009t iKuroda et al.ll2012t iTakiwaki et al.1 
12012ft . the enormous computational expense has lim- 
ited the most sophisticated supernova models, including 
the multi-species, multi-group neutrino transport, to 2D 
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axisymmetric simulations (|Ott et al.1 2008; Muller et al.1 
I2012T ) . Unfortunately, the fundamentally 3D hydrody- 
namics in the post-shock turbulent flow is qualitatively 
different if axisymmetry is imposed, as we discuss in de- 
tail below. Understanding the differences between 2D 
and 3D behavior and how dimension effects the mecha- 
nism of explosion is critical in the interpretation of re- 
alistic 2D simulations and ultimately in elucidating how 
real stars explode. If we find explosions in axisymmetry, 
should we expect them in 3D? If so, are the conditions 
identified in 2D as crucial to producing explosions man- 
ifest in 3D? How are explosions triggered in 3D? 

In this work, we do not perform sophisticated radia- 
tion hydrodynamic simulations of realistic core-collapse 
supernovae. Rather, we perform a series of simpli- 
fied numerical experiments designed to clarify how di- 
mension effects the hydrodynamics leading to explosions 
(see alsolMurphv fc Burrowsll2008t INordhaus et al.1120101: 
IHanke et al.l I2012T ) . We find that there are many differ- 
ences between 2D and 3D models, some quite dramatic, 
and that conditions identified as important in 2D may 
not remain so in more realistic 3D models. Nevertheless, 
our 3D models explode whenever the 2D models explode 
and, moreover, 3D models explode earlier. 

We begin in Sec. [5] with a discussion of our numerical 
setup and solution technique. Section[3]gives an overview 
of the basic results of our simulations and discusses qual- 
itative differences in the structures of 2D and 3D models. 
Section [4] begins the quantitative analyses of the simu- 
lations, showing how the global structures of the flows 
are different between 2D and 3D models. Section [5] com- 



2 



Dolence et al. 



pares various measures of the 2D and 3D turbulence in 
the post-shock flows. Section [6] discusses popular explo- 
sion metrics and introduces a simple model for explosions 
based on the runaway growth of bubbles. Finally, Sec. [7] 
discusses our results and conclusions. 

2. NUMERICAL SETUP 

Our setup is the sa me as that p resen ted in 

iBurrows et all (|2012fl and iMurphv et all (j2012t ). We 
use the CASTR O adaptive mesh r efinement hydro- 

" f20T0T) 



dynamics code (jAlmgren et al.l 120101) to evolve two- 
dimensional (axisymmetric) and three-dimensional mod- 
els of the collapse, bounce, and subsequent evolution 
of the 15- M^ non-rotating solar-metallicity progenitor 
of iWooslev &: Weaver! (|1995l ). The adaptive mesh uses 
six levels of factor two refinement with sa 0.5 km reso- 
lution in the inner 50 km and 2 km or better resolution 
everywhere behind and including the shock during the 
stalled pre-explosion phase. We ensure that both two- 
and three-dimensional simulations utilize the same re- 
finement criteria to minimize any differences that might 
arise from different grid structures. The domains include 
the inner 5000 km in radius. 

As in previous works, we adopt the "light bulb" 
prescription for ne utrino heating and cooling 
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(IMurphv fc Burrows! 120081: iNordhaus et all 
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2010; 



Hanke et al.l 12012! .IBurrows i et all l2012t IMurphv et all 
2012t iCouchl 12013 ). In this prescription, neutrino 



heating is parametrized by a constant driving electron 
neutrino luminosity L Ue (the electron neutrino and an- 
tineutrino luminosities are assumed to be equal) and we 
present results for three luminosities: 2.1 x 10 52 ergs _1 , 
2.2 x 10 52 ergs _1 , and 2.3 x lO^ergs" 1 . Neutrino 
cooling occurs at a rate oc T 6 (IBethdll990h . Since we do 
not explicitly treat the neutrino transport, the electron 
fra ction Y K is evolved accordi ng to the prescription given 
by iLiebendorfer et ail (|2005[ ) . Finally, the equation of 
state is ba s ed on the relativistic mean-field theory of 
IShen et al.1 (11998allbh and we assume nuclear statistical 
equilibrium. 

3. OVERVIEW OF SIMULATION RESULTS 

We find that the structure and evolution of the post- 
bounce hydrodynamics can be clearly distinguished be- 
tween simulations that differ only in dimension. In this 
section, we develop a qualitative view of the structure 
and evolution of 2D and 3D models. 

Figure Q] shows the development of the average shock 
radius for all three driving neutrino luminosities in both 
2D and 3D. For the first 100 ms after bounce, the mod- 
els are nearly indistinguishable, but diverge thereafter. 
The stalled shock radii in 3D are generally larger than in 
2D and are less variable in both angle and time. The 
models with driving luminosities of 2.2 x 10 52 ergs _1 
and 2.3 x 10 52 ergs _1 have sufficient neutrino heating 
to explode within the simulated time in both 2D and 
3D. The explosions, however, occ ur earlier in 3D than in 
2D, in qualitative agreement with lNordhaus et al.l (120101) 
despite deficiencies in that work ([Burrows et al.l [2012T ). 
Moreover, once explosions set in, the shock radii grow 
monotonically in 3D whereas they show oscillatory be- 
havior in 2D, at least early in the explosion phase. 

To develop a qualitative sense for some of the underly- 
ing hydrodynamic differences, we show snapshots of var- 
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Figure 1. Average shock radii for all six models considered in 
this work. The very early phases are quite similar, but the models 
diverge after 100 ms post-bounce. In the quasi-steady accretion 
phase, the average shock radius is ~30 km larger in 3D (solid) than 
in 2D (dashed) for a given driving neutrino luminosity. Finally, 
when explosions occur, they set in ~ 100-300 ms earlier in 3D than 
in 2D. 

ious quantities at select times in the evolution. Figure [5] 
shows slices of the entropy, radial velocity, magnitude of 
the vorticity (|V x v\), and velocity divergence (V • v) 
at 250 ms post-bounce for the L Ve — 2.2 x 10 52 ergs _1 
models in 2D and 3D. Figure [3] shows the same quan- 
tities at 500 ms post-bounce. Perhaps the most obvious 
difference between 2D and 3D is the clear presence of a 
preferred axis, leading to a distinctly prolate distortion of 
the post-shock flow. This is a generic result seen in all 2D 
simulations, even those that include so phisticated neu- 
trino transport and relativistic eff ects (lOtt et al.l 120081 : 
iMarek fc Jankal l2009t iMiiller et all |2012(). Importantly, 
this may be an artifact of assuming axisymmetry, the 
consequences of which are difficult to clarify. Another 
important difference, identified previously, is the exis- 
tence of more sma ll-scale structure in the flow in 3D 
(|Hanke et al.ll2012l) . as can be seen in the entropy, radial 
velocity, vorticity, and velocity divergence. Interestingly, 
when comparing the entropy and radial velocity maps 
in both the 2D and 3D models, there is a strong corre- 
spondence between large-scale structures; high-entropy 
plumes are associated with outflow whereas low-entropy 
regions are associated with inflow, a natural consequence 
of buoyancy-driven convection, which it has been argued 
dominates the flow in the stalled accretio n shock phase 
(|Burrows et alj|2012t IMurphv et al.ll2qi2D . In 3D, these 
rising plumes also have lower vorticities and larger ve- 
locity divergences than the surrounding flow, suggesting 
that these structures are relatively coherent and expand- 
ing. In 2D, these associations are more difficult to make 
by eye. 

The 2D and 3D post-shock flows have qualitative dif- 
ferences that are easily identified. The geometries of the 
flows are different and the characteristics of the turbu- 
lent, convective, post-shock flows are different. One con- 
sequence of these differences is a larger averaged stalled 
shock radius in 3D relative to 2D. In the remainder of this 
work, we undertake detailed analyses of the post-shock 
hydrodynamics to better understand how these qualita- 
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Figure 2. Snapshots of the entropy, radial velocity, magni- 
tude of the vorticity (|V X v\), and velocity divergence (V ■ v) 
at 250 ms post-bounce for the 2D (left) and 3D (slice, right) 
L„ e = 2.2 X 10 52 ergs -1 models. As evident in all four quantities 
and in contrast with the 2D models, 3D models show significantly 
more small scale structure, have no preferred axis, and tend to be 
more spherical in the early stalled accretion shock phase. 

tive differences manifest quantitatively in quantities that 
have been suggested to be of importance. 

4. GLOBAL STRUCTURES 

4.1. Integral Quantities 

While unable to capture the full scope of the differ- 
ences between 2D and 3D flows, integral quantities offer 
the advantages of simplicity and widespread usage. In 
the context of the neutrino mechanism, the integrated 
heating and cooling in the post-shock flow are naturally 
important quantities. Figures |H O and [6] show the total 




Figure 3. Same as Fig. [2] but at 500 ms post-bounce. As com- 
pared with the snapshot at 250 ms, the 3D model is developing 
significant asymmetry as it evolves towards explosion. The 2D 
model is qualitatively similar at 500 ms post-bounce to 250 ms post- 
bounce, with its distinctive prolate distortion and characteristically 
larger structures compared with the 3D model. 

integrated heating, cooling, and net heating minus cool- 
ing rates, respectively, in the region between the neu- 
trinosphere (where the optical depth to neutrinos is ap- 
proximately unity) and shock. While there tends to be 
marginally more heating in the 3D models, they have 
significantly more cooling and, therefore, less net heat- 
ing than their corresponding 2D models. If we focus on 
the gain region (i.e. only the region with net neutrino 
heating), we see in Fig. [7] that the 2D models have sig- 
nificantly more heating than their 3D counterparts. This 
arises, in part, because, prior to explosion, there is more 
mass in the gain region in 2D than in 3D, as shown in 



4 



Dolence et al. 



0.8 




0.0 0.2 0.4 0.6 0.8 1.0 

Time after bounce [s] 



Figure 4. Integrated rate of neutrino heating in the volume be- 
tween the neutrinosphere (r„) and shock (r s ). The rates are quite 
comparable between 2D (dashed) and 3D (solid), but 3D tends to 
be a few percent (~1— 3%) higher. 
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Figure 5. Integrated rate of neutrino cooling in the volume be- 
tween the neutrinosphere (rv) and shock (r B ). The 3D (solid) mod- 
els tend to cool more (~10%) than their 2D (dashed) counterparts. 

Fig. [HI Finally, in spite of the higher heating rates in 2D, 
the average specific entropy in the gain region, 

(s)gain = j^— I psdV, (1) 

-'"gain J gain 

is larger i n 3D than in 2D, as s hown in Fig. |9] and sug- 
gested bv lNordhaus et al~l (I2010D . This is consistent with 
the results shown in lrlanke et al.l ()2012h who found some- 
what larger average entropies in 3D than in 2D. We note, 
however, that a higher average entropy is not directly 
related to earlier explosions, as is clearly demonstrated 
by the 3D L Ue = 2.1 x 10 52 ergs -1 model, which has the 
highest average entropy of any model shown but does not 
explode within the simulated time. In any case, since en- 
tropy depends on the integrated heating, not the heating 
rate, this suggests that perhaps some material is exposed 
longer to heating in 3D than in 2D, to which we return 
to in § EUJ 
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Figure 6. Net rate of neutrino heating in the volume between the 
neutrinosphere (rv) and shock (r a ). Note that the heating rates are 
all negative, indicating that there is net neutrino cooling overall. In 
spite of the marginally higher heating rate, the 3D models (solid) 
show less (~10%) net heating per unit time than the 2D models 
(dashed) . 
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Figure 7. Net neutrino heating rate in the gain region between 
the gain radius (r 9 ) and the shock (r s ). After an initial transient 
phase lasting ~50ms, 2D models (dashed) show ~30% higher heat- 
ing rates in this region than the 3D models (solid) for a given driv- 
ing luminosity. The specific heating rates, found by dividing the 
heating rates shown above by the mass in the gain region, are also 
larger in 2D than in 3D. 

4.2. Radial Profiles 

Moving beyond integral quantities, we can look at ra- 
dial profiles, averaged over solid angle, to distinguish 
the global structures of 2D and 3D models. In or- 
der to highlight the trends with luminosity and dimen- 
sion, we compute two sets of profiles, one time-averaged 
from 200 ms to 300 ms post-bounce and another time- 
averaged from 450 ms to 550 ms post-bounce. The L„ c = 
2.3 x 10 52 ergs -1 models are excluded from the latter 
set of plots because they are already well into explosion. 
Time averaging is essential for the 2D profiles to mini- 
mize the large fluctuations that obscure their underlying 
quasi-steady structure. 
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Figure 8. Mass in the gain region for all six models in units of 
10 Mq. Prior to explosion, the gain mass is generally larger in 
2D (dashed) than in 3D (solid), which may contribute to the higher 
neutrino heating rates. 
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Figure 9. Average specific entropy in units of fc;,/baryon in the 
gain region. In the stalled shock phase, the average entropies are 
higher in 3D (solid) by about one unit compared to 2D (dashed). 
The turndown in the L„ e = 2.2 X 10 52 ergs _1 and L„ c = 2.3 X 
10 52 ergs -1 models is associated with the relatively lower entropy, 
unheated material that enters into the post-shock region in the 
exploding phase. 

Figure [TU] shows the time- and spherically-averaged en- 
tropy profiles. Consistent with Fig. [1] the entropy pro- 
files clearly show that the shock radii are systematically 
larger in 3D than in 2D prior to explosion. The 3D mod- 
els also have higher peak entropies than their correspond- 
ing 2D models, consistent with Fig. |H1 Interestingly, the 
peak entropies vary little over the luminosity range con- 
sidered for models with a given number of dimensions, 
but there is a clear distinction between the 2D and 3D 
models. On the other hand, the entropy profiles between 
—50 km and — 100 km are remarkably similar between all 
models. 

Figure [Til shows the profiles of radial velocity. As with 
the entropy profiles, the radial velocities in the post- 
shock region vary more with dimension than over the 



range of luminosities of models with a given number of 
dimensions. The magnitudes of the radial velocities in 
the post-shock flow are generally larger in 2D than in 3D. 
This remains true also when doing mass- weighted spher- 
ical averages, though the profiles are somewhat closer. 
The larger post-shock radial velocities in 2D are associ- 
ated with lower densities and, therefore, lower heating 
rates for radii beyond ~ 120-1 50 km as compared with 
3D. This has the effect of lowering the temperatures, 
which suppresses cooling and moves the gain radius in- 
ward (relative to 3D) towards higher densities and neu- 
trino fluxes. Figures [12] and [13] show the radial profiles 
of the net neutrino heating rate and bear out these ar- 
guments. In the end, the smaller gain radius in 2D leads 
to a larger integrated net heating rate in the gain re- 
gion in 2D as shown in Fig. [7] and a larger mass in the 
gain region as shown in Fig. [8] even though the densities 
and net heating rates are larger in 3D at radii beyond 
—120 km and — 150 km, respectively. 

Finally, Fig. [14] shows the radial profiles of the trans- 
verse kinetic energy. From 200 ms and 300 ms and be- 
tween — 100 km and the shock, the 2D models have nearly 
twice as much transverse kinetic energy as the 3D models. 
Again, the profiles vary much more with dimension than 
between the different luminosity models of the same di- 
mension. At this stage in the evolution, the turbulence in 
the post-shock flow is artificially vigorous in 2D axisym- 
metric models, consistent with multidimensional sim ula- 
tions of stellar convection (IMeakin fc ArnetH 120071 ) . At 
later times, the transverse kinetic energies in the 3D 
models have grown to be roughly comparable with the 
2D models, though the L v = 2.2 x 10 52 ergs _1 models 
still differ by -30%. 

5. TURBULENCE DIAGNOSTICS 
5.1. Dwell-Time Distributions 

During the stalled accretion shock phase, fluid ele- 
ments advect through the gain region and eventually set- 
tle onto the proto-neutron star. In spherical symmetry, 
the time to advect through the gain region (the dwell 
time) is short and shared by all fluid elements belonging 
to the same mass shell. In multiple dimensions, the as- 
pherical shock structure and post-shock turbulence lead 
to a distribution of dwell times, with some fraction of 
the mass being exposed lo nger to net neutrino heating 
(|Murphv fc Burrowsl[200l . This effect can increase the 
neutrino heating efficiency and has naturally been sug- 
gested t o be one key element of the multidimensiona l 
picture (jBurrows et al.lll995t fMurphv fc Burrowsll2008D . 
Here, we focus on the L Ve — 2.1 x 10 52 ergs -1 models 
since the dwell time distribution is difficult to interpret 
once fluid elements begin to participate in explosion. 

We measure the dwell-time distributions by follow- 
ing the trajectories of Lagrangian tracer particles which 
move according to the time-dependent velocity field. We 
initialize the particles at a radius of 400 km at 250 mt0 
post-bounce and begin integrating their dwell times once 
they pass through the shock. In 3D, we use approxi- 
mately 2 18 particles distributed quasi-uniformly on the 
sphere. In 2D, the fewer degrees of freedom means that 

3 We also looked at results based on particles injected at 500 ms 
post-bounce and the basic conclusions remain unchanged. 
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Figure 10. Time- and spherically-averaged entropy profiles in the inner 500 km for 2D (dashed) and 3D (solid) models. The left panel 
shows the profiles time-averaged from 200-300 ms post bounce while the right panel shows the profiles time-averaged from 450-500 ms 
post-bounce. The L„ e = 2.3 X 10 52 ergs _1 models explode before 500ms and are therefore not included in the latter. The 3D models 
generically have higher averaged entropies between ~ 100 km (near the gain radius) and the shock. 




Figure 11. The same as Fig. 1101 but for the time- and spherically-averaged radial velocity. The 3D models (solid) tend to have smaller 
radial velocities between ~ 100 km and the shock than the 2D models (dashed). 



there are fewer independent trajectories for fluid ele- 
ments of a given mass shell, even at the same resolution. 
To combat this we use 16 shells of approximately 2 12 par- 
ticles injected over 2 ms for a total of 2 16 particles. In 
2D, the particles are distributed uniformly in angle and 
we weight each particle's contribution to the dwell time 
distribution by sin# to account for the fact that each 
particle represents an angle-dependent annulus. 

The results are shown in Fig. [T5] We find that there 
are two important effects. First, the mean dwell time, 
(r), is longer in 2D. In steady state, the mass in the 
gain region is related to the mean dwell time simply by 
Mg a i n = M (t) and so a longer mean dwell time in 2D is 
consistent with the larger gain mass seen previously in 
Fig. [8] Second, the dwell time distribution in 3D has a 
longer, shallower tail than the 2D distribution. Figure IT6l 
shows one minus the cumulative distribution functions 



from which it can be seen that about 25% of the material 
spends more time in the gain region in 3D than in 2D. 
By recording the peak entropy reached by each tracer 
particle, we have found that there is a strong correlation 
between dwell time and peak entropy, suggesting that 
these long lived trajectories may be responsible, at least 
in part, for the larger average entropies seen in the 3D 
models and discussed previously. 

5.2. Turbulent Energy Spectra 

Perhaps the single most distinguishing characteristic 
between the 2D and 3D turbulent post-sho ck flows is 
found in their energy spectra. As shown by iKraichnanl 
( 1967) and later confirmed experimental ly in various con- 
texts (see, e.g., iBoffetta fc Eckel [20121) . turbulent cas- 
cades are different in 2D and 3D. In 3D turbulence, en- 
ergy is the only constant of the motion and this leads to a 
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Figure 12. The same as Fig. 1101 but for the time- and spherically averaged net heating rate. Overall, the 2D (dashed) and 3D (solid) 
models have very similar net heating profiles. Note that the "gain radius" is found at ~ 100 km where the net heating changes sign. 
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Figure 13. The same as Fig. 1121 but zoomed into the region of significant net neutrino heating. The 2D models (dashed) tend to have 
more heating at smaller radii (near the gain radius where the net heating changes sign), but less heating at larger radii, compared to the 
3D models (solid). 



single turbulent cascade that transfers energy from some 
driving wavenumber kd towards larger k (smaller scales) . 
In 2D, both energy and squared vorticity are constants 
of the motion which leads to cascades of energy and en- 
strophy (proportional to the mean squared vorticity). 
The enstrophy cascade transports enstrophy in fc-space 
from the driving wavenumber kd toward larger k (smaller 
scales). This leads to a characteristic k~ 3 scaling of th e 
velocity energy spectrum for k > kd (|Kr aichnanl 1 1 9671 ) . 
The energy cascade, by contrast, transports energy from 
kd to smaller k (larger scales) and leads to a fc~ 5 / 3 scaling 
of the velocity energy spectrum for k < kd, in direct anal- 
ogy with the Kolmogorov theory of turbulence. This is 
the inverse energy cascade of 2D turbulence, which tends 
to exaggerate motions on the largest scales of the flow. 
These turbulence theories were developed within highly 
idealized setups, assuming, for example, steady isotropic 
turbulence, and do not necessarily apply directly to the 



turbulence seen in t he core-collapse con text. Neverthe- 
less, we find, as did IHanke et al.l (|2012t ). that the basic 
predictions of these theories — the predominance of en- 
ergy at the largest scales in 2D, the excess energy at 
the largest scales in 2D relative to 3D, and the shallower 
slope of the velocity energy spectrum for k > kd (and 
therefore more energy for large k) in 3D relative to 2D — 
are all confirmed in our simulations. 

The most natural basis to represent the matter fields 
in the quasi-spherical post-shock flow is the basis of 
real spherical harmonics. We decompose the arbitrary 
scalar quantity Q into spherical harmonics with time- 
and radially-dependent coefficients 



ai n (t, r) = i Q(t, r, 9, 0)y, m (0, 4>)dSi, (2) 
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Figure 14. The same as Fig. 1101 but for the time- and spherically-averaged transverse kinetic energy per unit mass. In the early stalled 
accretion shock phase (200-300 ms post-bounce), the 2D models (dashed) have ~50— 100% more specific transverse kinetic energy than the 
3D models (solid). This early turbulent vigor seems to be an artifact of assuming axisymmetry. At later times, the 2D/3D difference is 
reduced, but is still significant (~30%). 
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Figure 15. Distribution of dwell times in the gain region for tracer 
particles injected at 250 ms post-bounce in the 2D (blue) and 3D 
(red) L v<1 = 2.1 X 10 52 ergs -1 models. The 2D model has a longer 
mean dwell time, but the 3D distribution has a long shallow tail. 
Note that the data beyond 400 ms for the 2D model is dominated 
by shot noise associated with the finite number of tracer particles. 



where 



y/2 N™P{ n (cos 6)cosm<j) m > 0, 



Yr(e,4>) = { W (cos 0) 



m = 0, (3) 
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In 2D, axisymmetry implies that all coefficients with 
m =/= are identically zero. Here, we consider Q — 
{p,P,s,y/~pVi\, where Vi represents the spherical veloc- 
ity components {v r , vg} in 2D and {v r , vg, v^} in 3D. We 
compute the discrete energy spectrum as a function of 
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Figure 16. One minus the cumulative distribution function of 
dwell times in the gain region for 2D (blue) and 3D (red). The 
curves cross at ~50 ms, with the 3D values larger at later times, 
indicating that ~25% of the accreted material spends more time 
in the gain region in 3D compared to 2D. 



spherical harmonic degree I as 



(5) 



m— — I 



where it should be understood th at ai m , and, therefor e, 
E(l), depend on time and radius (jBurrows et al.ll20i"2l ). 

Figure [17] shows the energy spectra of yfpve for the 
L„ c = 2.1 x 10 52 ergs _1 models at a radius of 150km 
and time-averaged between 450-500 ms after bo unce. We 
also reproduce the results of lHanke et al.l (|2012l ) for com- 
parison. While these curves change in detail for dif- 
ferent quantities, over time, and at different radii, the 
qualitative trends and relation between 2D and 3D are 
quite robust. At I — 1, 2D has an order of magnitude 
more energy in all quantities Q — {p, P, s, ^/pvi}, consis- 



Dimensionality and the Hydrodynamics of CCSNe 



9 




10-' 



2D y/pVQ 
3D y/pve 



2D Hankc et al. (2012) 
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Figure 17. Discrete energy spectra of ^fpvg as a function of spherical harmonic degree I for the L„ E = 2.1 X 10 52 ergs 1 2D (blue) and 
3D (red) models, measured at r = 150 km and time-averaged between 450-500 ms post-bounce. The left panel includes power-laws <x Z -2,6 
and <x i -1 , indicating the inertial range scaling in 2D and 3D, respectively. The spectra show that the 2D mode l has excess power for 
all I < 30. Grid-scale dissipation sets in at I ~ 40. The right panel reproduces the results shown in Hankc et al. (2012) for comparison 
(dashed). Both set s of results show the same inertial range scalings, but the excess power at low I is even more extreme in the results from 
IHanke et aLl l2012T) . 



tent with the idea of an inverse energy cascade which 
pumps energy into the largest scales. The results of 



IHanke et alj ((2012) are even more extreme, with a fac- 
tor « 50 more energy in the I — 1 mode in 2D relative 
to 3D. In real-space, this excess energy at the largest 
scales manifests as the characteristic "sloshing" always 
found in 2D simulations (including those with coupled 
neutrino transport) . This sloshing i s typically associated 
with the development of the SASI dBlondin et al.l 120031: 
iScheck et alll2008t iFoglizzo et al.l 120071 ). but the inverse 
energy cascade will always produce excess energy at I = 1 
in 2D, even if the turbulence is driven by, for example, 
neutrino-driven convection. The SASI may be capable, 
in principle, of producing significant energy in I = 1, 
but disentangling its effects from the inevitable I = 1 
energy associated with the inverse cascade would seem 
to require 3D simulations. Whether the low-Z energy is a 
result of the SASI or not, all of the 3D supernova simula- 
tions thus far present ed in the literature show muted or 
nonexistent sloshing (llwakami et al.l 120081 : IHanke et al.l 
l2012t iBurrows et al.l I2012I ). though there has yet to be 
a fully self-consistent 3D radiation hydrodynamic simu- 
lation. Therefore, the results seen thus far suggest that 
these violent sloshing motions are an artifact of assum- 
ing axisymmetry, not a feature t hat must be incorpo rated 
into 3D models as suggested by IHanke et al.l (|2012f) . Im- 
portantly, this artifact is not a small effect; most of the 
energy in the flow in 2D simulations is at low-Z. This 
artifact, in fact, dominates the flow and has poorly un- 
derstood consequences on other coupled aspects of the 
problem, including the neutrino transport. 

At intermediate I, the spectra are significantly steeper 
in 2D (~ l~ 2 - 6 ) than in 3D (~ l' 1 ). That these power- 
laws differ from those naively expected from simple the- 
oretical arguments should not be surprising given that 
the turbulence analyzed here is, among other things, not 
steady-state or isotropic. Qualitatively, however, the ex- 
pectation of a steeper slope in 2D relative to 3D is con- 



firmed a nd we find our resu lts generally consistent with 
those of IHanke et alj (|2012j ffl The transition to these 
slopes occurs around I ~ 10 in 2D and I ~ 4 in 3D and re- 
flects a characteristic scale for convective plumes, a scale 
that likely depends on the model-dependent size of the 
gain region. The 3D model has significantly more energy 
at small-scales than the 2D model. At I > 40 (spatial 
scales ~10km at a radius of 150km), we begin to see the 
e ffects of grid- s cale (2 km at this radius) dissipation. 

IHanke et al.l (|2012l ) showed that their results, espe- 
cially in 3D, were sensitive to resolution. The comparison 
of energy spectra in Fig. [17] affords us the opportunity 
to directly compare the effective resolutions of our inde- 
pendent calculations by comparing the scales at which 
dissipation begins to set in. As noted above, our 3D re- 
sults begin to devi ate from the I -1 pow er-law at I s» 40, 
while the results in IHanke et al.l (|2012l ) begin to deviate 
at I w 25, confirming our expectation tha t our models 
have n early double the effective resolution. Han ke et al.1 
argue that 3D models become less prone to ex- 
plosion as resolution is increased, yet our higher resolu- 
tion 3D models explode earlier than the 2D models. The 
source of this apparent discrepancy is unclear. 

6. MULTIDIMENSIONAL EXPLOSIONS 

6.1. Explosion Conditions 

A number of quantities have been proposed in the lit- 
erature that are meant to distinguish exploding from 
non-exploding models and, further, to define in a sys- 
tematic way the time at which a model transitions into 
the exploding phase. The most widely used condition 
is based on the idea of a critical ratio of the advection 
to heating time scales dJanka fc KeilHl998t iThompsonl 
[2000HThompsonet all 120031) . There" are numerous ways 

4 IHanke et"aT] $20T3) report scalings of l~ 3 and l~ 5 / 3 in 2D and 
3D, respectively, but Fig. H7l shows that their results are consistent 
with the shallower slopes reported here. 
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of defining these timescales. Here, we define the advec- 
tion time as 

Jr. («r) 

where (v r ) is the spherically-averaged radial velocity. In 
multi-D models, the shock radius R s and the gain radius 
-Rgain are not uniquely defined boundfl We appeal to 
the radial profiles shown previously and define the shock 
radius as the outermost zero in the radial velocity gra- 
dient and the gain radius as the first zero crossing in 
the net heating rate interior to the shock (always around 
^100 km). Furthermore, to minimize the large fluctu- 
ations in t a( jv that appear in 2D due to transient and 
localized fluctuations in v r , we time- average the velocity 
profiles over a ±15 ms window. We define the heating 
timescale as 

_ SZJ^ M e °«P>> (Ye)))^rHr 

^hcat — ~7? 1 CO 

/£. ( P )((n-c))^dr 

•* "-gam 

where angle brackets indicate solid-angle averaging, e is 
the specific internal energy, e is f ne zero-point energy 
of the EOS, and T~L — C is the net heating rate per unit 
mass. The results are shown in Fig. [T5] After an initial 
transient phase, the models settle on an approximately 
(model-dependent) constant value. Rather than there 
being a particular critical value for the ratio of advection 
to heating time scales, explosions seem to be robustly 
connected to rapid growth from these model-dependent 
quasi-steady values. Nonetheless, we define an explosion 
time in a systematic manner by measuring the time at 
which this ratio exceeds 0.5 without later dropping below 
this value. 




Time after bounce [s] 



Figure 18. Ratio of advection to heating timescales (denned in 
the text) for the 3D (solid) and 2D (dashed) models. Explosions are 
associated with a rapid growth in this ratio, though it is difficult 
to identify a particular critical value for the onset of explosions. 
With a critical value of 0.5, the L Vt , = 2.3 X 10 52 ergs — 1 model 
explodes 159 ms earlier in 3D relative to 2D, while the L Vc = 2.2 X 
10 52 ergs -1 model explodes 333 ms earlier in 3D. 



5 Indeed, the gain region need not even be bounded from below 
by a single closed surface. 




Time after bounce [s] 

Figure 19. Maximum of the ratio of the square of the local 
sound speed to escape speed for the 3D (sol id) and 2D (dashed) 
model s, smoothed for clarity, as discussed by Pcicha & Thompson 
(2012) in their antcsonic condition. In agreement with lMiiller et all 
(2012), we find th a t the critical value of 0.2 suggested by 
Pcicha & Thompson (2012) is too low, but using a larger value 
(~ 0.3) may be a useful diagnostic of explosion. With a critical 
value of 0.3, the L Ua = 2.3 X 10 52 ergs -1 model explodes 87ms 
earlier in 3D relative to 2D, while the L„ c = 2.2 X 10 52 ergs _1 
model explodes 309 ms earlier in 3D. 

An alternative explosion condition (the "antesonic " 
condition) was suggested bv lPeicha fe Thompson! (|2012f) . 
Based on parametrized ID steady-state models, they sug- 
gested that explosions occur when max(c^/«j r ) r eaches 
a critical value w 0.2. IMurphv fc Meakinl (|2011[ ) 
tested this condition with their parametrized 2D sim- 
ulations and found it to be consistent with their results. 
IMuller et~a l. (2012) on the other hand, using 2D radia- 
tion hydrodynamic simulations, assessed the validity of 
the antesonic condition and suggest that it may not be a 
robust indicator of explosion, but, in any case, the crit- 
ical value should at lea st be somewhat larger (~0.35). 
Like IMuller et~aTl ((2012T ). we find that a larger value for 
the critical condition is required and we adopt 0.3 as 
the critical ratio. Unfortunately, uiax(c^/v^ sc ) is a noisy 
quantity, with brief spikes as large as sa 0.5 that then re- 
turn well below the critical value. The smoothed curves, 
however, seem to reliably indicate when explosions set 
in, but these smoothed curves are necessarily produced 
ex post facto. In other words, the antesonic condition is 
not a reliable indicator of explosion given the instanta- 
neous value of the ratio. The evolutions for the 2D and 
3D models, smoothed for clarity, are shown in Fig. fT9l 
As above, we identify the time of explosion as the time 
when the smoothed version of max(c 2 /-i> 2 sc ) exceeds 0.3 
without later dropping below this value. 

A final, phcnomcnological, explosion condition is sim- 
ply when the average shock radius exceeds 400 km 
without later recedin g below this value, as used in 
iNordhaus et al.l (|2010f ). In this case, explosion times can 
be read directly from the shock radius evolution curves 
in Fig. [Q 

Table [1] shows the explosion times and correspond- 
ing mass accretion rates as determined by the three 
conditions above for the L Ve = 2.2 x 10 52 ergs _1 and 
L Uc = 2.3 x 10 52 ergs -1 models. All three conditions give 
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Table 1 

Explosion Times and Accretion Rates 



r a 

L "e 


Dimension 






Explosion Condition 










^adv /^heat 


max 


(csA4c) 




(Rs) 






(ms) 


M 

(MqS" 1 ) 


^cxp 

(ms) 


M 

(Mes- 1 ) 


^cxp 

(ms) 


M 

(Mgs- 1 ) 


2.2 


2D 
3D 


1031 
698 


0.182 
0.213 


920 
611 


0.196 
0.224 


954 
695 


0.188 
0.214 


2.3 


2D 
3D 


571 
412 


0.231 
0.269 


501 
414 


0.241 
0.268 


557 
434 


0.233 
0.262 



a 10 52 ergs _1 

comparable numbers, but the antesonic condition tends 
to give the earliest indication of explosion. We note, how- 
ever, that this conclusion may be somewhat sensitive to 
the particular critical values adopted. That the three 
conditions give comparable explosion times is not sur- 
prising; the ratio of advection to heating times, the ratio 
of sound to escape speed, and the average shock radius all 
grow as models transition into explosion. None of these 
conditions is able to robustly predict when or even if an 
explosion will occur before the explosion begins. This 
suggests that while these conditions may be indicative of 
explosion, they are by no means the full story. 

With the results shown in Table [T] in hand, we can 
quantify the delay between the 2D explosions and the 
earlier 3D explosions. Taking the minimum difference 
found between the explosion conditions, we find that the 
L Vll = 2.3 x 10 52 ergs _1 3D model explodes at least 87 ms 
before the corresponding 2D model, while the L„ e = 
2.2 x 10 52 ergs _1 3D model explodes at least 259 ms ear- 
lier than the corresponding 2D model. These delays may 
translate into ap preciable differences in the explosion en- 
ergy at infinity (jYamamoto et al1l2Q12f ). with earlier ex- 
plosions (3D) being more energetic. 

6.2. Explosion Trigger in 3D 

What conditions trigger explosions in 3D? The multi- 
D nature of the hydrodynamics is generally agreed to 
be important in producing such conditions. Turbu- 
lence in the post-shock flow leads to a distribution of 
dwell times for accreting parcels of matter, increasing 
the mass in the gain region and, therefore, the efficiency 
of neutrino heating. In this way, turbulence plays a cru- 
cial, but secondary , role as an aid to neutrino heating. 
IMurnhv fc Meakinl j20Tl|) argue that convection modi- 
fies the quasi-steady global structure of the flow by in- 
troducing turbulent fluxes of, for example, enthalpy and 
entropy. Their model focuses on how turbulent convec- 
tion modifies averaged radial profiles, rather than on the 
effects of particular turbulent fluctuations, and was able 
to account for the differences in, for example, the radial 
profiles of entropy between ID and 2D models. While 
these are important roles for the post-shock turbulent 
motions, here we suggest that the turbulent fluctuations 
themselves, driven by neutrino heating, may be instru- 
mental in triggering explosions. 

The dominant hydrodynamic instability, aside from 
the explosion itself, in our param etrized 3D models is 
neutrino-heating-dr iven convection (|Burrows "elall[2012l: 
iMurphv et al.ll2012D . In the nonlinear phase, the post- 
shock turbulent flow involves a complex interaction be- 



tween buoyantly rising, neutrino-heating-driven plumes, 
negatively buoyant accretion streams, turbulent entrain- 
ment resulting from secondary Kelvin-Helmholtz in- 
stabilities, and dissipative interactions between rising 
plumes and the bounding shock. The evolution of a 
buoyant "bubble" depends on parameters, including the 
heating rate, bubble size, and background inflow velocity. 
Small bubbles tend to be shredded by turbulent entrain- 
ment, while large bubbles can rise all the way to the shock 
and locally push it outward. Our simulations suggest 
that some bubbles can continue rising, locally pushing 
the shock to ever larger radii until the global explosion 
is triggered. Figure [20] shows a sequence of volume ren- 
derings illustrating that large buoyant features form in 
the flow and push the shock to larger radii. Figure [21] 
shows the evolution of the shock surface where the same 
local growth of the shock radius can be clearly seen. In- 
terestingly, the growth of these features occurs not on a 
dynamical time, but appears to proceed quasi-statically. 
The same basic picture holds in all of our 3D simula- 
tions, suggesting that this may be a generic feature of 
the transition to explosion. 

We can try to understand the conditions for the run- 
away growth of bubbles (and the associated triggered ex- 
plosions) by considering the highly simplified toy model 
presented in the Appendix. In this model, the evolu- 
tion of a bubble's outer radius (Rb) is determined by 
the competition between neutrino driving power (L v t, 
L v is the neutrino driving luminosity and r is the ef- 
fective optical depth of the bubble) and accretion power 
(aGMM I Rb, a is a constant defined in the Appendix) 
associated with the ram pressure of material immediately 
behind the shock. The evolution follows 



a 







dR b 

~~dT 4tt GMM b 



L,,t — a 



gmm\ 

R b I 



which has the solution 



R b {t) 



aGMM 



where 



L v t + e xt (aGMM/R - L v t) ' 
fl M 



A = a- 



' 4tt M b ' 



(8) 



(9) 



(10) 



Here, Slo is the bubble's constant solid angle and Mb is 
the bubble's fixed mass. Qualitatively, R b (t) has two 
types of behavior. When the neutrino driving term ex- 
ceeds the ram pressure term, Rb{t) grows without bound, 
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Figure 20. Sequence of volume renderings of the specific entropy for the 3D L„ e = 2.2 X 10 ergs - model. As seen at the bottom left 
of each image, a high entropy plume forms around ~250-300ms and persists for hundreds of milliseconds, eventually pushing the shock 
out far enough to seemingly trigger the global explosion. A similar structure appears around ~450 ms at the top right of each image, which 
leads to similar local shock expansion thereafter. 

Second, there is a critical luminosity for a given mass 
accretion rate and shock radius above which a bubble 
will runaway. If the stalled shock radius scales as L@, 
then the critical luminosity for runaway bubble growth 
is proportional to M 1 / 1+ ' 3 . Empirically, j3 ~ 3 and if we 
adopt the parameters used above we find 



presu mably triggering the global explosion ()Thompsonl 
2000). When the ram pressure term dominates, the bub- 
ble recedes, perhaps allowing another bubble to take its 
place. 

This simple model makes two predictions. First, the 
model predicts a characteristic growth of the maximum 
shock radius, at least in the quasi-static growth phase. 
Figure 1221 shows the maximum shock radii (i? max ) from 
the L Ue = 2.2 x KFergs" 1 and L Ue = 2.3 x lO^ergs" 1 
3D models along with fitted model solutions. To com- 
pute the fits, we fix M = 1.6 M , M = 0.25 MqS" 1 , 
and Mb/fio = 2.5 x 10 30 g and fit the model solution to 
the data between 250 ms post-bounce and the time when 
-Rmax = 1000km, leaving Rq, t, and a as free param- 
eters. Both fits give r ps 0.04 and a ps 0.25. Given 
the simplicity of the model, the agreement between the 
hydrodynamical simulation results and the model predic- 
tions is quite surprising and encouraging. 



L 



crit 



i2.2 



M 



1/4 



M 



1/4 



1.6 M 

T \-l/4 



.0.04 



\ 0.25 M Q s- 
x 10 52 ergs -1 , 



(11) 



which seems consistent with the critical explosion curves 
from parametrized multi-D models shown in other 
works (iMurphv fc Burrows! 120081 ; [ Nordhaus et~al1 l2010t 
IHanke et al.ll2012t ICouchl 120121) . Inasmuch as runway 
bubble growth is associated with explosions, this may 
be viewed as an alternative, albeit crude, derivation of 
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Figure 21. Evolution of the shock surface from the 3D L„ e = 2.2 X 
aspherical nature of the shock surface and the secularly growing feature 



the critical luminosity curve of lBurrows fc Goshvl {1993). 
This may suggest that the reduction in the critical lumi- 
nosity in going from ID to 2D and 3D models might 
arise, in part, from the emergence of bubbles and their 
runaway growth. 

7. DISCUSSION AND CONCLUSIONS 

We have presented analyses of parametrized 2D ax- 
isymmetric and 3D core-collapse supernova models. Our 
basic conclusion is not surprising nor controversial — the 
hydrodynamics of core-collapse supernovae are different 
between 2D axisymmetric and 3D models. These dif- 
ferences are many and, while some are subtle and per- 
haps not crucial to understanding the mechanism, some 
are quite dramatic and make interpreting 2D supernova 
models problematic. 

Our parametrized models indicate that the global 
structures of the flows are different between 2D and 3D. 



10 erg s model. The colors indicate radius to em pha size the 
associated with the high-entropy plume shown in Fig. 1201 

We see this reflected in integral measures like the mass 
and average entropy in the gain region and in the radial 
profiles of various quantities. In spite of identical heat- 
ing and cooling prescriptions between 2D and 3D, our 
analyses show how the different global structures effect 
the heating and cooling rates. We find that the 2D mod- 
els have significantly higher integrated net heating rates 
than their corresponding 3D models and attribute this 
to the smaller gain radius (and more mass in the gain 
region) in 2D, which may in turn be related to the larger 
radial velocities in 2D. On the other hand, the 3D mod- 
els tend to have higher densities at large radii and larger 
average shock radii. 

The dwell time distributions for accreted parcels of 
matter offer another vantage point from which to distin- 
guish 2D and 3D models. In both cases, the post-shock 
turbulence leads to a broad distribution of dwell times. 
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Figure 22. Maximum shock radii for the 3D Z/„ e = 2.2 X 
10 52 ergs -1 and L„ c = 2.3 X 10 52 ergs -1 models, along with fitted 
solutions (black dashed lines) based on Eq. [9] See the text and 
Appendix for discussions. 

In our simulations, the mean dwell time is somewhat 
longer in 2D than in 3D, consistent with the larger mass 
in the gain region in 2D than in 3D, but the 3D distri- 
bution has a relatively prominent tail towards long dwell 
times. About 25% of the material spends more time in 
the gain region in 3D than in 2D, being exposed to more 
integrated heating and reaching higher peak entropies. 

Ultimately, many of the differences we see are plausibly 
associated with the character of the post-shock turbulent 
flow. In 2D, turbulent energy is pumped into the largest 
scales of the flow, which inevitably gives rise to the slosh- 
ing behavior manifest in all modern 2D core-collapse su- 
pernova simulations. This sloshing behavior has been 
identified as a crucial ingredient in nearly all successful 
neutrino-driven explosion models to date (iBuras et al.1 
I20M IMarek fc Jankal IMiiller et al.H2012h . In 3D~ 
this sloshing motion is muted or absent, at least in part, 
because the turbulent energy transport is predominantly 
towards small scales. Since we see explosions earlier in 
3D than in 2D, vigorous sloshing is either not critical 
in any dimension or the explosion mechanism operates 
differently in 2D and 3D. 

Finally, we present a toy model that describes the 
evolution of buoyant bubbles driven by the competition 
between neutrino heating power and accretion power. 



The simple model is able to account for the quasi- 
static growth of the maximum shock radius in the ex- 
ploding 3D models with neutrino luminosities of L„ c = 
2.2 x 10 52 ergs _1 and L Ve = 2.3 x lO^ergs -1 . It also 
predicts the existence of a critical luminosity for a given 
mass accretion rate and shock radius beyond which a 
bubble will have runaway growth. We speculate that this 
runaway growth triggers the global explosion of the star, 
and therefore that the critical luminosities for runaway 
bubble growth and explosion are effectively the same. 

The mean background and fluctuating turbulent com- 
ponents of the post-shock flow in 2D axisymmetric mod- 
els are different from 3D models and, more importantly, 
likely not representative of the hydrodynamics of super- 
nova cores in Nature. While the models presented in this 
work are incomplete, lacking adequate neutrino trans- 
port and feedback, this basic conclusion seems robust. 
It seems inevitable that we must await 3D models be- 
fore drawing conclusions concerning the trigger of core- 
collapse supernovae. 
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APPENDIX 

TOY MODEL OF A NEUTRINO-DRIVEN BUBBLE 

Consider a bubble of fixed solid angle fio pinned beneath the shock at radius Rb- The energy injection rate due 
to neutrino heating is ^(r2n/47r)L„T. For a quasi-steady system, this energy injection rate should be balanced by 
the rate of energy loss. The interaction between the bubble and the surrounding flow will remove energy from 
the bubble at a rate ~ Cdpv 3 floRl, where Cd is an effective drag coefficient that encapsulates information about 
the geometry and character of the flow, including effects like turbulent entrainment. For a bubble pinned beneath 
the shock, the density and velocity are given by the shock jump conditions, so this loss term can be rewritten as ~ 
C d {Q /ATT)(f 2 /n 2 )(2GMM/R b ), where / is the pre-shock fraction of the free-fall velocity (w 1/V2) and K is the shock 
compression ratio. Also, lifting the bubble out of the gravitational potential requires a power ~ dRb/dtGMMb/R^. 
Equating the driving neutrino power to the sum of the ram and gravitational powers, we find 

dR b n Q R 2 b ( gmm\ 

= — L u t — a Al 

dt 4-7T GMMb \ R b J y ' 



Dimensionality and the Hydrodynamics of CCSNc 
where a = 2Cdf 2 fR- 2 is a dimensionless constant of order unity. 
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